%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Computing steady state equilibrium with persistent idiosyncratic
%%%% productivity

clc; clear; close all; format short


%% 1. ------------------------ Control of experirements--------------------
%%%% control
comparative_statics      = 1;
comparative_statics_plot = 1;

THETA_grid          = [0.45, 0.50, 0.55];
THETA_grid_legend   = {'0.45', '0.50', '0.55'};

options = optimset('Display','on'); 
options.OptimalityTolerance = 1e-5; 
options.StepTolerance = 1e-5;
options.MaxIter = 50000; 
options.MaxFunEvals = 50000; 
options.FunctionTolerance = 1e-5;
options.TolX = 1e-5; 
options.TolFun = 1e-5;

colors              = [0, 0, 1; 1, 0, 0; 0.9, 0.6, 0.2];
marker_style        = {'none'; 'none'; 'none'};
line_style          = {'--';'-';'-.'};

%%%% choose either gap_a or persistence
variable_x          = 'gap';
%variable_x         = 'persistence' 
gap_a_grid          = linspace(0, 0.0234, 5); % this controls whether the persistent case is the same as the IID case. Set it to zero to have the IID result
prob_grid           = [0.5, 0.65, 0.75]; % this controls whether the persistent case is the same as the IID case. Set it to zero to have the IID result
prob_1              = 0.75;  % transition matrix persistence; prob of staying in high state
prob_2              = 0.75;


%% 2.1 ----------   Solve cases with idiosyncratic shocks   ---------------

%%%% parameters
load param_calibrated % from calibration when THETA = 0.50;
par                     = par_calibrated;
par_iid                 = par;
par_persist             = par;
par_persist.endo_entry  = 0; % no endogenous entry

disp('*********************************************************************')
disp('I am solving IID cases to be used later')


x0       = [par_iid.EPS_L_BOUND * 4 + 1, 0.15, 0.11];
initial_ = zeros(7, length(THETA_grid));
for ii = 1: length(THETA_grid)
    par_iid.THETA       = THETA_grid(ii);   
    fun                 = define_function_fun(par_iid); 
    x                   = fsolve(@(x)eqm_ss_iid_case_6(x, par_iid, fun, 1, 0), x0);   
    [fval, eqm]         = eqm_ss_iid_case_6(x, par_iid, fun, 1, 0);    
    L              = x(1);
    B              = x(2);
    ALPHA_0        = x(3);
    initial_(:,ii) = [eqm.Z; eqm.Z; eqm.K; eqm.K; eqm.ALPHA_0; eqm.ALPHA_0; eqm.B];
end


%% 2.2 -----------     Do persistent-case computation  --------------------
disp(' ')
disp('*********************************************************************')
fprintf('I am solving the case of "%s" for persistent productivity shocks \n',... 
variable_x)


switch variable_x
    case 'gap'
        l_x = length(gap_a_grid);
        x_grid = gap_a_grid;
    case 'persistence'
        l_x = length(prob_grid); 
        x_grid = prob_grid;
end
l  = length(THETA_grid);
Z_1_result              = zeros(l, l_x);
Z_2_result              = zeros(l, l_x);
Z_result                = zeros(l, l_x);
K_1_result              = zeros(l, l_x);
K_2_result              = zeros(l, l_x);
K_result                = zeros(l, l_x);
ALPHA_0_1_result        = zeros(l, l_x);
ALPHA_0_2_result        = zeros(l, l_x);
ALPHA_0_result          = zeros(l, l_x);
C_result                = zeros(l, l_x);
Inv_result              = zeros(l, l_x);
H_result                = zeros(l, l_x);
Y_result                = zeros(l, l_x);
Prod_result             = zeros(l, l_x);
R_ratio_result          = zeros(l, l_x);
P_share_result          = zeros(l, l_x);
Prob_target_1_result    = zeros(l, l_x);
Prob_target_2_result    = zeros(l, l_x);
Prob_aqc_1_result       = zeros(l, l_x);
Prob_aqc_2_result       = zeros(l, l_x);
Prob_aqc_11_result      = zeros(l, l_x);
Prob_aqc_12_result      = zeros(l, l_x);
Prob_aqc_21_result      = zeros(l, l_x);
Prob_aqc_22_result      = zeros(l, l_x);
Prob_aqc_bb_result      = zeros(l, l_x);
Prob_aqc_bs_result      = zeros(l, l_x);
Prob_aqc_sb_result      = zeros(l, l_x);
Prob_aqc_ss_result      = zeros(l, l_x);
aqc_1_result            = zeros(1, l_x);
aqc_2_result            = zeros(1, l_x);
HT_credit_1_result      = zeros(1, l_x);
HT_credit_2_result      = zeros(1, l_x);


%%%% Define indicators
%%%% 1 stands for high state; 2 stands for low state; 1 meets 1, 1 meets 2, 2 meets 1, and 2 meets 2
%%%% the 2nd column means the current state; the 1st column means the previous state

par_persist.state = [repmat([1,1], 4, 1); repmat([1,2], 4, 1); repmat([2,1], 4, 1); repmat([2,2], 4, 1)]; 
par_persist.state_tilde     = repmat([1, 1; 1, 2; 2, 1; 2, 2], 4, 1);

position            = 2 * (par_persist.state(:,1) - 1) + par_persist.state(:,2); %locate the position in P for eps type
position_tilde      = 2 * (par_persist.state_tilde(:,1) - 1) + par_persist.state_tilde(:,2); % locate the position in P for eps_tilde type 

par_persist.H_yes   = par_persist.state(:,1) == 1; % extract indicator for high state yesterday;
par_persist.L_yes   = par_persist.state(:,1) == 2; % extract indicator for low state yesterday;
par_persist.H_tod   = par_persist.state(:,2) == 1; % extract indicator for high state yesterday;
par_persist.L_tod   = par_persist.state(:,2) == 2; % extract indicator for low state yesterday;
par_persist.yesterday = [par_persist.H_yes'; par_persist.L_yes'];
par_persist.today   = [par_persist.H_tod'; par_persist.L_tod'];

%%%% conditional on yesterday's status or today's status
%%%% selecting different possibilities of pairs
select_ = par_persist.yesterday;
select_ = [select_; 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;...
    0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0;...
    0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0;...
    0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1];

bb = (par_persist.state(:,1)==1).* (par_persist.state_tilde(:,1)==1); % big buy big
bs = (par_persist.state(:,1)==1).* (par_persist.state_tilde(:,1)==2); % big buy small
sb = (par_persist.state(:,1)==2).* (par_persist.state_tilde(:,1)==1); % small buy big
ss = (par_persist.state(:,1)==2).* (par_persist.state_tilde(:,1)==2); % small buy small;

tic
if comparative_statics == 1
    
    for ii = 1: length(THETA_grid)
        fprintf('I am doing computation for theta = %.2f', THETA_grid(ii))
        %THETA_grid(ii)
        par_outer           = par_persist;
        par_outer.THETA     = THETA_grid(ii);
        fun                 = define_function_fun(par_outer);
        Z_1_temp            = zeros(1, l_x);
        Z_2_temp            = Z_1_temp;
        K_1_temp            = Z_1_temp;
        K_2_temp            = Z_1_temp;
        ALPHA_0_1_temp      = Z_1_temp;
        ALPHA_0_2_temp      = Z_1_temp;
        Prob_target_1_temp  = Z_1_temp;
        Prob_target_2_temp  = Z_1_temp;
        Prob_aqc_1_temp     = Z_1_temp;
        Prob_aqc_2_temp     = Z_1_temp;
        Prob_aqc_11_temp    = Z_1_temp;
        Prob_aqc_12_temp    = Z_1_temp;
        Prob_aqc_21_temp    = Z_1_temp;
        Prob_aqc_22_temp    = Z_1_temp;
        Prob_aqc_bb_temp    = Z_1_temp;
        Prob_aqc_bs_temp    = Z_1_temp;
        Prob_aqc_sb_temp    = Z_1_temp;
        Prob_aqc_ss_temp    = Z_1_temp;
        aqc_1_temp          = Z_1_temp;
        aqc_2_temp          = Z_1_temp; 
        HT_credit_1_temp    = Z_1_temp;
        HT_credit_2_temp    = Z_1_temp;
        Z_temp              = Z_1_temp;
        K_temp              = Z_1_temp;
        Y_temp              = Z_1_temp;
        C_temp              = Z_1_temp;
        Inv_temp            = Z_1_temp; 
        H_temp              = Z_1_temp; 
        Prod_temp           = Z_1_temp;
        R_ratio_temp        = Z_1_temp;
        P_share_temp        = Z_1_temp; 
        ratio = 1.1;
        x0                  = initial_(:,ii); 
        x0(1) = x0(1) / ratio; x0(2) = x0(2) * ratio; 
        x0(3) = x0(3) * ratio; x0(4) = x0(4) / ratio;

        for jj = 1: l_x

            a = zeros(2,1); 
            P = zeros(2,2); 
            par_temp = par_outer; 

            switch variable_x
                case 'gap'
                    a(1)    = log(1 + gap_a_grid(jj)); 
                    a(2)    = log(1 - gap_a_grid(jj));
                    P(1,1)  = prob_1; 
                    P(1,2)  = 1 - prob_1; 
                    P(2,1)  = 1 - prob_2; 
                    P(2,2)  = prob_2;
                case 'persistence'
                    a(1)    = log(1 + gap_a);
                    a(2)    = log(1 - gap_a);
                    P(1,1)  = prob_grid(jj); 
                    P(1,2)  = 1 - prob_grid(jj); 
                    P(2,1)  = 1 - prob_grid(jj); 
                    P(2,2)  = prob_grid(jj);
            end

            eigenvector              = null(eye(size(P,1)) - P');
            par_temp.P               = P;  
            par_temp.Pss             = eigenvector / sum(eigenvector);     
            par_temp.a               = a;
            par_temp.q_average_prod  = arrayfun(@(a_)integral(@(x) x .* fun.f(x, a_), fun.l_bound(a_), fun.h_bound(a_)), a); % the average productivity        
            par_temp.Pss_meet        = par_temp.Pss(par_temp.state(:,1));
            par_temp.Pss_tilde_meet  = par_temp.Pss(par_temp.state_tilde(:,1));
            P_trans                  = P';
            par_temp.P_tr_meet       = P_trans(position);
            par_temp.P_tr_tilde_meet = P_trans(position_tilde);

            x = fsolve(@(x)eqm_ss_persistent_case_6_log_normal(x, par_temp, fun, 0), x0, options);

            [fval, c_star, K, Y, G, Z, H, K_bar, E_a, E_p, q_capital, q_money, q_profit, w, ...
                L_s, L_b, ALPHA_0, ALPHA_0_meet, ALPHA_0_tilde_meet, Pss_tilde_meet, K_meet, K_tilde_meet, a_meet, a_tilde_meet,...
                int_aqc, int_sppe, int_prob_aqc, int_prob_target, int_average_2nd_price, int_prod_std, int_HT_credit, int_debt] =...
                eqm_ss_persistent_case_6_log_normal(x, par_temp, fun, 1);

            fval;
            Z_1_temp(jj)          = Z(1);
            Z_2_temp(jj)          = Z(2);       
            K_1_temp(jj)          = K(1);
            K_2_temp(jj)          = K(2);
            ALPHA_0_1_temp(jj)    = x(5);
            ALPHA_0_2_temp(jj)    = x(6);

            Prob_aqc_temp           = select_ * ( par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_aqc);
            Prob_target_temp        = select_ * (ALPHA_0_meet .* par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_target);
            aqc_temp                = select_ * (ALPHA_0_meet .* par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* K_meet .* int_aqc);   
            HT_credit_temp          = select_ * (ALPHA_0_meet .* par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_HT_credit); 

            Prob_aqc_1_temp(jj)     = Prob_aqc_temp(1); 
            Prob_aqc_2_temp(jj)     = Prob_aqc_temp(2);
            Prob_aqc_11_temp(jj)    = Prob_aqc_temp(3)/par_temp.P(1,1);
            Prob_aqc_12_temp(jj)    = Prob_aqc_temp(4)/par_temp.P(1,2);
            Prob_aqc_21_temp(jj)    = Prob_aqc_temp(5)/par_temp.P(2,1);
            Prob_aqc_22_temp(jj)    = Prob_aqc_temp(6)/par_temp.P(2,2);

            Prob_target_1_temp(jj)  = Prob_target_temp(1); 
            Prob_target_2_temp(jj)  = Prob_target_temp(2);

            Prob_aqc_bb_temp(jj)    = bb' * (par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_aqc);
            Prob_aqc_bs_temp(jj)    = bs' * (par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_aqc);
            Prob_aqc_sb_temp(jj)    = sb' * (par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_aqc);
            Prob_aqc_ss_temp(jj)    = ss' * (par_temp.P_tr_meet .* par_temp.P_tr_tilde_meet .* Pss_tilde_meet .* int_prob_aqc);

            aqc_1_temp(jj)          = aqc_temp(1);
            aqc_2_temp(jj)          = aqc_temp(2);
            HT_credit_1_temp(jj)    = HT_credit_temp(1);
            HT_credit_2_temp(jj)    = HT_credit_temp(2);


            Z_temp(jj)              = par_temp.Pss(1) * Z(1) + par_temp.Pss(2) * Z(2);
            K_temp(jj)              = par_temp.Pss(1) * K(1) + par_temp.Pss(2) * K(2);
            Y_temp(jj)              = Y;
            C_temp(jj)              = c_star;
            Inv_temp(jj)            = par_temp.Pss' * K * par_temp.DELTA;
            H_temp(jj)              = H;
            Prod_temp(jj)           = K_bar^par_temp.ETA;
            R_ratio_temp(jj)        = (E_a + E_p) / (E_a + E_p + Inv_temp(jj));
            P_share_temp(jj)        = E_p / (E_a + E_p);
            
            x0 = x; %%%% used as the next initial value; which means that we can parfor this loop        
        end

        Z_1_result(ii,:)            = Z_1_temp;
        Z_2_result(ii,:)            = Z_2_temp;
        K_1_result(ii,:)            = K_1_temp;
        K_2_result(ii,:)            = K_2_temp;
        ALPHA_0_1_result(ii,:)      = ALPHA_0_1_temp;
        ALPHA_0_2_result(ii,:)      = ALPHA_0_2_temp;

        Prob_target_1_result(ii,:)  = Prob_target_1_temp;
        Prob_target_2_result(ii,:)  = Prob_target_2_temp;
        Prob_aqc_1_result(ii,:)     = Prob_aqc_1_temp;
        Prob_aqc_2_result(ii,:)     = Prob_aqc_2_temp;
        Prob_aqc_11_result(ii,:)    = Prob_aqc_11_temp;
        Prob_aqc_12_result(ii,:)    = Prob_aqc_12_temp;
        Prob_aqc_21_result(ii,:)    = Prob_aqc_21_temp;
        Prob_aqc_22_result(ii,:)    = Prob_aqc_22_temp;
        Prob_aqc_bb_result(ii,:)    = Prob_aqc_bb_temp;
        Prob_aqc_bs_result(ii,:)    = Prob_aqc_bs_temp;
        Prob_aqc_sb_result(ii,:)    = Prob_aqc_sb_temp;
        Prob_aqc_ss_result(ii,:)    = Prob_aqc_ss_temp;

        aqc_1_result(ii,:)          = aqc_1_temp;
        aqc_2_result(ii,:)          = aqc_2_temp;   

        HT_credit_1_result(ii,:)    = HT_credit_1_temp;
        HT_credit_2_result(ii,:)    = HT_credit_2_temp;
        Z_result(ii,:)              = Z_temp;
        K_result(ii,:)              = K_temp;
        Y_result(ii,:)              = Y_temp;
        C_result(ii,:)              = C_temp;
        Inv_result(ii,:)            = Inv_temp;
        H_result(ii,:)              = H_temp;
        Prod_result(ii,:)           = Prod_temp;
        R_ratio_result(ii,:)        = R_ratio_temp;
        P_share_result(ii,:)        = P_share_temp; 

    end

end
toc


%% 3. plot comparative statics
 
if comparative_statics_plot == 1
    x_grid = x_grid * 2;    
    xq = linspace(0, x_grid(end), length(x_grid) * 2);
    
    int_method = 'spline';
    
    for ii = 1:l
        Yq(ii,:)        = interp1(x_grid, Y_result(ii,:), xq, int_method);
        Iq(ii,:)        = interp1(x_grid, Inv_result(ii,:), xq, int_method);
        Cq(ii,:)        = interp1(x_grid, C_result(ii,:), xq, int_method);
        Hq(ii,:)        = interp1(x_grid, H_result(ii,:), xq, int_method); 
        R_ratio_q(ii,:) = interp1(x_grid, R_ratio_result(ii,:), xq, int_method); 
        P_share_q(ii,:) = interp1(x_grid, P_share_result(ii,:), xq, int_method); 
        Z_1_q(ii,:)     = interp1(x_grid,Z_1_result(ii,:), xq, int_method); 
        Z_2_q(ii,:)     = interp1(x_grid,Z_2_result(ii,:), xq, int_method);
        K_1_q(ii,:)     = interp1(x_grid,K_1_result(ii,:), xq, int_method); 
        K_2_q(ii,:)     = interp1(x_grid,K_2_result(ii,:), xq, int_method); 
        Prob_target_1_q(ii,:) = interp1(x_grid,Prob_target_1_result(ii,:), xq, int_method); 
        Prob_target_2_q(ii,:) = interp1(x_grid,Prob_target_2_result(ii,:), xq, int_method);
        Prob_aqc_1_q(ii,:)    = interp1(x_grid,Prob_aqc_1_result(ii,:), xq, int_method); 
        Prob_aqc_2_q(ii,:)    = interp1(x_grid,Prob_aqc_2_result(ii,:), xq, int_method);  
        
        Prob_aqc_bb_q(ii,:)   = interp1(x_grid,Prob_aqc_bb_result(ii,:), xq, int_method); 
        Prob_aqc_bs_q(ii,:)   = interp1(x_grid,Prob_aqc_bs_result(ii,:), xq, int_method); 
        Prob_aqc_sb_q(ii,:)   = interp1(x_grid,Prob_aqc_sb_result(ii,:), xq, int_method); 
        Prob_aqc_ss_q(ii,:)   = interp1(x_grid,Prob_aqc_ss_result(ii,:), xq, int_method); 
        
        
        liq_output_h_q(ii,:)  = interp1(x_grid, (par.CHI_k * K_1_result(ii,:) + ALPHA_0_1_result(ii,:) .* Z_1_result(ii,:) + HT_credit_1_result(ii,:)) ./ Y_result(ii,:) * 100, xq, int_method); 
        liq_output_l_q(ii,:)  = interp1(x_grid, (par.CHI_k * K_2_result(ii,:) + ALPHA_0_2_result(ii,:) .* Z_2_result(ii,:) + HT_credit_2_result(ii,:)) ./ Y_result(ii,:) * 100, xq, int_method); 
        cash_output_h_q(ii,:) = interp1(x_grid, ALPHA_0_1_result(ii,:) .* Z_1_result(ii,:) ./ Y_result(ii,:) * 100, xq, int_method); 
        cash_output_l_q(ii,:) = interp1(x_grid, ALPHA_0_2_result(ii,:) .* Z_2_result(ii,:) ./ Y_result(ii,:) * 100, xq, int_method); 
    end

    figure()
    
    subplot(3,4,1)
    h = plot(xq * 100, K_1_q, 'Linewidth', 2); grid on
    title('Capital (H group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,3])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);
    
    subplot(3,4,2)
    h = plot(xq * 100, K_2_q, 'Linewidth', 2); grid on
    title('Capital (L group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,3])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);

    subplot(3,4,3)
    h = plot(xq * 100, cash_output_h_q, 'Linewidth', 2); grid on
    title('Cash/output (H group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylabel('%')
    ylim([3,5])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);
    
    subplot(3,4,4)
    h = plot(xq * 100, cash_output_l_q, 'Linewidth', 2); grid on
    title('Cash/output (L group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylabel('%')
    ylim([3,5])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);
    
    subplot(3,4,5)
    h = plot(xq * 100, liq_output_h_q, 'Linewidth', 2); grid on
    title('Liquidity/output (H group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylabel('%')
    ylim([3,5])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);
   
    subplot(3,4,6)
    h = plot(xq * 100, liq_output_l_q, 'Linewidth', 2); grid on
    title('Liquidity/output (L group)','Fontsize',14)
    xlabel(append('%',variable_x))
    ylabel('%')
    ylim([3,5])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);
      
    subplot(3,4,7)
    h = plot(xq * 100, R_ratio_q, 'Linewidth', 2); grid on
    title('R share','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0.15,0.55])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);

    subplot(3,4,8)
    h = plot(xq * 100, P_share_q, 'Linewidth', 2); grid on
    title('P share','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0.15,0.55])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);

    subplot(3,4,9)
    h = plot(xq * 100, Prob_aqc_bb_q, 'Linewidth', 2); grid on
    title('Prob to acquire: BB','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,0.25])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);

    subplot(3,4,10)
    h = plot(xq * 100, Prob_aqc_bs_q, 'Linewidth', 2); grid on
    title('Prob to acquire: BS','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,0.25])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);    
    
    subplot(3,4,11)
    h = plot(xq * 100, Prob_aqc_sb_q, 'Linewidth', 2); grid on
    title('Prob to acquire: SB','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,0.25])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);

    subplot(3,4,12)
    h = plot(xq * 100, Prob_aqc_ss_q, 'Linewidth', 2); grid on
    title('Prob to acquire: SS','Fontsize',14)
    xlabel(append('%',variable_x))
    ylim([0,0.25])
    set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style, {'Marker'}, marker_style, 'MarkerSize', 5);      
     
    for i = 1:length(THETA_grid)
        %entry_legend{i} =  ['\theta = ' num2str(THETA_grid(i))];
        entry_legend{i} =  ['\theta = ' THETA_grid_legend{i}];
    end
    legend(entry_legend, 'location','Southeast','Orientation','horizontal','Fontsize',14);
    
 
end

